function y = noise_energy (chi,span,CONSTS)
% if chi==1 - chi(t) is the Gaussian white noise time dependence
% if chi==0 - chi(t) is a pulsed time dependence

    a = CONSTS.a;
    Vn = CONSTS.Vn;
    c = CONSTS.c;
    Z0 = 4*pi/c;
    d = CONSTS.d;
    
    wLH = CONSTS.omega_LH;
    wH  = CONSTS.omega_H;
    omega1 = span(1); %wLH; %%1.21e-2*wH;
    omega2 = span(2); %wH-5; %1e-4*wH;
    
    delta_w= omega2-omega1;
    num_points = size(chi,1);
    omega_vec = linspace(wLH,wH,num_points)'; %(omega1:(omega2-omega1)/100:omega2);
    chi_from_w = func_chi_from_omega(omega_vec,CONSTS,false);
    if(false)
        func_plot = func_Wmn(omega_vec,chi_from_w,1,1,CONSTS);
        figure; plot(omega_vec,func_plot);
    end
    
    n_vec = (1:300)';
    m_vec = (0:20)';    
    
    matlabpool;
  
    for mi = 1:size(m_vec, 1)
        m = m_vec(mi);
        
        parfor ni = 1:size(n_vec, 1)
            n = n_vec(ni);

            if(chi)
                Wmn(mi,ni) = 8*(pi/delta_w)^2*((Vn^2)/(Z0^2*(log(4*a/d))^2))*quadgk(@(omega)(((func_chi_mean_from_omega(omega',CONSTS,false))').^2 ...
                                                         .*func_Wmn_without_chi(omega,m,n,CONSTS)),omega1,omega2);
            else
                Wmn(mi,ni) = 8*(pi/delta_w)^2*((Vn^2)/(Z0^2*(log(4*a/d))^2))*quadgk(@(omega)(...
                                                         func_Wmn(omega,1,m,n,CONSTS)),omega1,omega2);
            end
        end  
        
        Wm(mi) = sum(Wmn(mi,:));    
        fprintf('m= %d   Wm= %e\n', m,Wm(mi));
    end       
    
    matlabpool close;
    y = 2*sum(Wm);
end

